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The spontaneous emergence of coherent behavior through synchronization plays a key role in neural 
function, and its anomalies often lie at the basis of pathologies. Here we employ a parsimonious 
(mesoscopic) approach to study analytically and computationally the synchronization (Kuramoto) 
dynamics on the actual human-brain connectome network. We elucidate the existence of a so-far-uncovered 
intermediate phase, placed between the standard synchronous and asynchronous phases, i.e. between order 
and disorder. This novel phase stems from the hierarchical modular organization of the connectome. Where 
one would expect a hierarchical synchronization process, we show that the interplay between structural 
bottlenecks and quenched intrinsic frequency heterogeneities at many different scales, gives rise to 
frustrated synchronization, metastabUity, and chimera-like states, resulting in a very rich and complex 
phenomenology. We uncover the origin of the dynamic freezing behind these features by using spectral 
graph theory and discuss how the emerging complex synchronization patterns relate to the need for the 
brain to access -in a robust though flexible way- a large variety of functional attractors and dynamical 
repertoires without ad hoc fine-tuning to a critical point. 



N euro-imaging techniques have allowed the reconstruction of structural human brain networks, com- 
posed of hundreds of neural regions and thousands of white-matter fiber interconnections. The 
resulting "human connectome" (HC)''^ turns out to be organized in moduli -characterized by a 
much larger intra than inter connectivity- structured in a hierarchical nested fashion across many scales^^ '°. 
On the other hand, "functional" connections between nodes in these networks have been empirically inferred 
from correlations in neural activity as detected in electroencephalogram and functional magnetic resonance 
time series. Unveiling how structural and functional networks influence and constrain each other is a task of 
outmost importance. A few pioneering works found that the hierarchical-modular organization of structural 
brain networks has profound implications for neural dynamics'* "^'''. For example, neural activity propagates 
in hierarchical networks in a rather distinctive way, not observed on simpler networks'^; beside the usual two 
phases -percolating and non-percolating- commonly encountered in models of activity propagation, an 
intermediate "Griffiths phase""" emerges on the hierarchical HC network'''" '". Such a Griffiths phase stems 
from the existence of highly diverse relatively-isolated moduli or "rare regions" where neural activity remains 
mostly localized generating slow dynamics and very large responses to perturbations'"" 

Brain function requires coordinated or coherent neural activity at a wide range of scales, thus, neural synchron- 
ization is a major theme in neuroscience^'"'^". Synchronization plays a key role in vision^', memory^^, neural 
communication^'', and other cognitive functions"^''. An excess of synchrony results in pathologies such as epilepsy or 
Parkinsonian disease, while neurological deficit of synchronization has been related to autism and schizophrenia^l 
Our aim here is to scrutinize the special features of synchronization dynamics"'" -as exemplified by the 
canonical Kuramoto modeP'"^'- running on top of the best available human connectome mapping''^ ™. This 
consists of a network of 998 nodes, each of them representing a mesoscopic population of neurons -able to 
produce self- sustained oscillations"- whose mutual connections are encoded by a symmetric weighted connec- 
tivity matrix W''^. The validity of this admittedly simplistic Kuramoto model as a convenient tool to explore the 
generic features of complex brain dynamics at a large scales has been recently emphasized in the literature^'"^^. 
Here, we uncover the existence of a novel intermediate phase for synchronization dynamics -similar in spirit to 
the Griffiths phases discussed above- which stems from the hierarchical modular organization of the HC and 
which gives rise to very complex and rich synchronization dynamical patterns. We identify this novel phase as the 
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Figure 1 | Global synchronization dynamics in the human connectome. (a) Time average of the order parameter R{t), for Kuramoto dynamics on the 
HC network for a specific and frxed set of frequencies extracted from a N{Q, 1) Gaussian distribution. A broad intermediate regime separates the 
incoherent phase (low k) from the synchronous one (high k). In this regime, coherence increases with k in an intermittent fashion, and with strong 
dependence on the frequency realization, (b) Raster plot of individual phases (vertical axis) showing local rather than global synchrony and illustrating the 
coexistence of coherent and incoherent nodes (k = 2.7) as time runs, (c) R{t) for 4 values of k (arrows in the main plot), (d) Adjacency matrix of the HC 
network with nodes ordered to emphasize its modular structure as highlighted by a community detection algorithm (main text), keeping the partition 
into the 2 hemispheres (dashed lines). Intra-modular connections (shown in color) are dense while inter-modular ones (grey) are limited to tiny subsets, 
acting as interfaces between moduli. Integration between hemispheres is mostly carried out by the 3 central moduli. This plot visually illustrates the 
hierarchical modular organization of the human connectome network. 



optimal regime for the brain to harbor complex behavior, large 
dynamical repertoires, and optimal trade-offs between local segrega- 
tion and global integration. 

The Kuramoto dynamics on a generic network (see™ for a nice and 
comprehensive review) is defined by: 

N 

m =(n, + kJ2 sin [0,(0 - 0,(0] , (1) 

where 0,(f) is the phase at node i at time t. The intrinsic frequencies oj, 
- accounting for node heterogeneity- are extracted from some arbit- 
rary distribution function g{a}), Wjj are the elements of the N X N 
connectivity matrix W, and k is the coupling strength. Time delays, 
noise, and phase frustration could also be straightforwardly imple- 
mented. The Kuramoto order parameter is defined as Z{t) = 
R{t)e"^^'^ = ^e'"'''' ^, where 0 £ R{t) £ 1 gauges the overall coherence 
and (//(f) is the average global phase. In large populations of well- 
connected oscillators without frequency dispersion, perfect coher- 
ence {R = I) emerges for any coupling strength; on the other hand, 
frequency heterogeneity leads to a phase transition at some critical 
value of k, separating a coherent steady state from an incoherent 
one^^"^". Analytical insight onto this phase transition can be obtained 
using the celebrated Ott-Antonsen (OA) ansatz, allowing for a pro- 
jection of the high-dimensional dynamics into an evolution equation 
for Z{t) with remarkable accuracy in the large-Nlimit^''"'\ 

Results 

Novel phase between order and disorder in the HC. We have 
performed a computational study of the Kuramoto model running 
on top of the HC network (details are given in the Methods section). 
Our results reveal the existence of an intermediate regime placed 
between the coherent and the incoherent phase (see Fig. 1). This is 
characterized by broad quasi-periodic temporal oscillations of R(t) 



which wildly depend upon the realization of intrinsic frequencies^''". 
Anomalously large sampling times would be required to extract good 
statistics for the actual mean values and variances. Collective oscil- 
lations of R{t) are a straightforward manifestation of partial syn- 
chronization and they are robust against changes in the frequency 
distribution (e.g. Gaussian, Lorentzian, uniform, etc.) whereas the 
location and width of the intermediate phase depend upon details. As 
this phenomenology is reminiscent of Griffiths phases -posed in 
between order and disorder and stemmig from the existence of 
semi-isolated regions'"" it is natural to investigate how the HC 
hierarchical modular structure affects synchronization dynamics. 

Any network with perfectly isolated and independently synchro- 
nized moduli trivially exhibits oscillations of R{t), with amplitude 
peaking at times when maximal mutual synchronization happens to 
be incidentally achieved. Such oscillations can become chaotic if a 
finite and relatively small number of different coherent moduli are 
coupled together ''*. Thus, in a connected network without delays or 
other additional ingredients, oscillations in the global coherence are 
the trademark of strong modular structure with weakly intercon- 
nected moduli. 

Strong modular organization into distinct hierarchical levels is 
indeed present in the HC as reveled by standard community detec- 
tion algorithms'^ "" and as already discussed in the literature (see e.g.'^ 
and references therein). For instance, we have found that the optimal 
partition into disjoint communities -i.e. the partition maximizing 
the modularity parameter'"'- corresponds to a division in 12 com- 
munities (see Fig. Id) while, at a higher hierarchical level, a separa- 
tion into just 2 moduli -the 2 cerebral hemispheres- is obtained^ 
(Fig. Id). Obviously these 2 coarser moduli include the 12 above as 
sub-moduli. Although more levels of hierarchical partitioning could 
be inferred (see e.g.'" and refs. therein), for the sake of simplicity we 
focus on these two levels I, I = 1 and I = 2 with 12 and 2 moduli, 
respectively. 
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Figure 2 | Two-block model, (a) Sketch of the two-block model, (b) Global order parameter for the two-block model with M = 128 and two interfacial 
nodes. Results of the numerical integration of the 258 Kuramoto equations (blue points) are in strikingly good agreement with the integration of 
Eqs.(2) (solid blue line). Local block-wise order parameters are shown for comparison (small symbols; dashed lines are guides to the eye). A first 
transition, where local order emerges, occurs at fc = 0.02, while global coherence is reached atk~ 90. In the intermediate region, R(f) oscillates (inset), 
revealing the lack of global coherence. Despite the simplicity of this toy model, these results constitute the essential building-block upon which further 
levels of complexity rely (see main text). 



A simplistic model for global oscillations. To shed further light on 
the properties of synchronization on the HC, we consider a very 
simple network model -allowing for analytical understanding- 
which will constitute the elementary "buUding-block" for subse- 
quent more complex analyses. This consists of a few blocks with 
very large internal connectivity and very sparse inter-connectivity. 
Each block is composed by a bulk of M^l nodes that share no 
connection with the outside and a relatively small "interfacial" set 
that connects with nodes in other blocks. For instance, in the simplest 
realization, consisting of just two blocks connected by a single pair of 
nodes (Fig. 2), each block is endowed with local coherence Ta.b, 
average phase i/za.B) and average characteristic frequency 0)^3, 
while 1-node interfaces have perfect coherence r = 1, phase (^a.b- 
and characteristic frequency \'a,b. In this case, N = 2M + 2, and the 
OA ansatz can be safely applied to each block (large M) but not to 
single-node interfaces. In the particular case (convenient for 

w) are zero-mean Lorentz distri- 



analytical treatment) in which 

butions e(oj) = s 

of OA equations can be easily shown to be: 



with spreads 8a.b> the resulting set 



l/'A = «A+fc 



l + ri 
2rA 



sm{VA-i^A 



-[MrA + cos(«!>A-i/'A)] 



(2) 



^A = Va + k[MrA sm{>pA -^Ia) +^H<Ps-<Pa)] 



(together with r = 1 for each 1-node interface), and a symmetric set 
(A B) for block B. The solution of Eq.(2) -displayed in Fig. 2- 
reveals a transition to local coherence within each block at a certain 
threshold value of A: = 0.02. As soon as local order is attained, Ta.b ~ 1 
and \I/as~^> from Eq.(2) the mutual synchronization process obeys 



9A''{^'A + MwA)+ksm{(pj^-(Pf^) 



(3) 



and a symmetrical equation for q>^. For small k, the right-hand side is 
dominated by Va + Moja' whereas the average value cua becomes ar- 
bitrarily small within blocks (assuming that M is large), the frequency 
Va does not. Consequently, synchronization between the two blocks 
through the interfacial link is frustrated: each block remains inter- 



nally synchronized but is unable to achieve coherence with the other 
over a broad interval of coupling strengths. This interval is delimited 
above by a second transition at fc ~ max{M|coA,B|) \'a,b}> where k is 
large enough as to overcome frustration and generate global 
coherence. This picture is confirmed by numerical integration of 
the full system of N coupled Kuramoto equations as well as by its 
OA approximation (Eq.(2)), both in remarkably good agreement. 
Therefore, local and global coherences have their onsets at two 
well-separated transition points^^ and -similarly to the much more 
complex HC case- _R oscillates in the intermediate regime (Fig. 2). 
Similar results hold for versions of the model with more than two 
moduli (e.g. 4; see below). The existence of two distinct (local and 
global) transitions had already been reported in a recent study of 
many blocks with much stronger inter-moduli connections than 
here'^ (even if, owing to this difference, no sign of an intermediate 
oscillatory phase was reported). In particular, the value of two-block 
models has already been explored in the past, for systems of identical 
oscillators with non-zero phase lags, in which each node is coupled 
equally to all the others in its community, and less strongly to those in 
the other^'. In such systems, local coherence emerged for large 
enough values of the phase lag. Our two-block model shows that 
the presence of "structural bottlenecks" between moduli combined 
with heterogeneous frequencies at their contact nodes (interfaces) 
are essential ingredients to generate a broad region of global oscilla- 
tions in R, even in the absence of phase lag. Still, it is obviously a too- 
simplistic model to account for all the rich phenomenology emerging 
on the HC, as we show now. 

Oscillations of local coherence in the HC. Fig. 3 shows numerical 
results for the local order parameter r'" for some of the moduli at the 
2 hierarchical levels, / = 1 and Z = 2 in the HC network. It reveals that 
(Fig. 3a) local coherences exhibit oscillatory patterns in time (with 
characteristic frequencies typically between 0.01 and 0.1 Hz) and 
that (Fig. 3b) the transition to local coherence at progressively 
higher hierarchical level occurs at progressively larger values of k; 
i.e. coherence emerges out of a hierarchical bottom-up process as 
illustrated above for the for the two-block model (see'^'*^). Observe, 
however, that local oscillations were not present in the two-block 
model. This suggests that the 12 moduli in the HC are on their turn 
composed of finer sub-moduli and that structural frustration, as 
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Figure 3 | Local synchronization in the human connectome. (a) Oscillations of the local order parameters ("chimera-like states") in one particular 
modulus in the partitions oftheHC into 12 (green, / = 1, and fc = 3) and 2 (magenta, / = 2, and k= W) moduli, respectively. The characteristic frequency 
of these oscillations is typically between 0.01 and 0.1 Hz (a range which coincides with slow modes detected in brain activity; seee.g.-'^). (b) Average of the 
local order parameter over all moduli and ( c) chimera index for moduli at levels as in a) , as a function of k. Global order (thin black line in b) ) emerges only 
after local order is attained at lower levels, (d) Average decay of activity p for identical frequencies co = 0 in the HC network and comparison with a single- 
level modular network (made up of 4 similar random moduli at a single hierarchical level) of the same size and average connectivity as the HC network. 
Symbols stand for different values of A:, (e) Characteristic decay times corresponding to the inverse of the first 1 000 non-trivial eigenvalues of the Laplacian 
matrix (xaxis) as a function of their respective ordered indices (y axis), for networks as in (d). The stretched exponential behavior in (d) is the result of the 
convolution of slow time scales associated with small eigenvalues in (e). 



introduced above, affects all hierarchical levels. The average variance 
of local coherences (called chimera index, xY* exhibits a marked peak 
-reflecting maximal configurational variability- at the transition 
point for the corresponding level (Fig. 3b-c and Methods section). 
Similar intra-modular oscillatory patterns -dubbed chimera states- 
have been recently found'" '''' ''^ in Kuramoto models in which explicit 
phase lags induce a different kind of frustration, hindering global 
synchronization. Strictly speaking, chimeras are defined in systems 
of identical oscillators. In such a case, a non-zero phase lag term is 
essential for partial synchronization to occur. Realistic models of the 
brain, however, require oscillators to be heterogenous. States of 
partial synchronization in empirical brain networks with frequency 
heterogeneity have been found for Kuramoto models with explicit 
time delays^'. In contrast, the chimera-like states put forward here 
have a purely structural origin, as they arise from the network 
topology. It was noted in the past that synchronization in a 
synthetic network with hubs could be limited to those hubs by 
tuning clustering properties, and global order could be attained in 
a monotonous step-like fashion upon increasing k*^. Fig. 3b instead 
reveals that the ordering process in the hierarchical modular HC may 
be non-monotonous: coherence does not systematically grow with k. 
Indeed, the emergence of local order in some community may hinder 
or reduce coherence in others, inducing local "desynchronization" 
and reflecting the metastable nature of the explored states. 

Anomalous dynamics on the HC. Fixing all intrinsic frequencies to 
be identical allows us to focus specifically on structural effects. Thus, 
we consider, without loss of generality, the simple case cu, = 0, and 
define the "activity" p = 1 — (R). In this case, perfect asymptotic 
coherence should emerge for all values of k but, as illustrated in 
Fig. 3d, the convergence towards p = 0 turns out to be extremely 
slow (much slower than exponential). This effect can be analytically 
investigated assuming that, for large enough times, all phase diffe- 



rences are relatively small. Then, up to first order, 0, = — . LyOj 
where L,, = Sjj'EiWji — are the elements of the Laplacian 
matrbc"''"". Solving the linear problem, = e"*^"v|v|0;(O), 

where denotes the /-th Laplacian eigenvalue (0 = Ai < I2 < . . . < 
Xf^ and v| the !-th component of the corresponding eigenvector. 
Given that the averaged order parameter can be written as 

2^)^ I averaging over initial conditions, 

and considering that (as the Laplacian has zero row-sums'"*) 2i = 
0, we obtain 

p(0 = yf:e-"", (4) 

^ ( = 2 

where o is the standard deviation of the initial phases. This 
expression holds for any connected network. As usual, the larger 
the spectral gap ~Ai, the more "entangled"''" the network and thus 
the more difficult to divide it into well separated moduli (I2 = 0 
only for disconnected networks)'"'**. For large spectral gaps all 
timescales are fast, and the last expression can be approximated by 
its leading contribution, ensuing exponential relaxation to p = 0, as 
in fact observed in weU-connected network architectures (Erdos- 
Renyi, scale free, etc.''"). This is not the case for the HC matrix, for 
which a tail of small non- degenerate eigenvalues is encountered (see 
Fig. 3e and'^). Each eigenvalue a, in the tail corresponds to a natural 
division of moduli into submoduli'"*, and the broad tail reflects the 
heterogeneity in the resulting modular sizes. As a consequence, each 
of these eigenvalues -with its associated large timescale, f, = 1/1 — 
contributes to the sum above, giving rise to a convolution of 
relaxation processes, entailing anomalously-slow dynamics, which 
could not be explained by a single-level modular network (see 
Fig. 3d-e): slow dynamics necessarily stems from the existence of a 



SCIENTIFIC REPORTS | 4:5990 | DOI: 1 0. 1 038/srep05990 



4 



hierarchy of moduU and structural bottlenecks. As explained in 
Methods, in the case of the HC the convolution of different times 
scales gives rise to stretched-exponential decay, which perfectly fits 
with numerical results in Fig. 3. It was noted in the past that strongly 
modular networks exhibit isolated eigenvalues in the lower edge of 
the laplacian spectrum. Synchronization would develop in a step- 
wise process in time, where each transient would be given by each 
isolated eigenvalue''"'. In our case, the depth of the hierarchical 
organization and the strength of topological disorder produce 
instead a quasi- continuous tail of eigenvalues, and the step-wise 
process is replaced by an anomalous stretched-exponential behavior. 

A more refined model: hierarchical modular synthetic networks. 

To shed additional light on the previous findings for the HC -i.e. the 
emergence of chimera-like states and anomalously slow dynamics- 
we suggest to go beyond the single-level modular network model and 
study hierarchical modular networks (HMN) in which moduli exists 
within moduli in a nested way at various scales' HMN are 
assembled in a bottom-up fashion: local fuUy- connected moduli 
(e.g. of 16 nodes) are used as building blocks. They are recursively 
grouped by establishing additional inter-moduli links in a level- 
dependent way as sketched in Fig. 4(top)'^''''. 

Our computational analyses of the Kuramoto dynamics on HMN 
substrates (see Fig. 4) reveal: (i) a sequence of synchronization tran- 
sitions for progressively higher hierarchical levels at increasing 
values of k, (ii) chimera-like states at every hierarchical level, result- 
ing in a hierarchy of metastable states with maximal variability at the 
corresponding transition points, (iii) extremely slow relaxation 
toward the coherent state when all internal frequencies are identical. 
Furthermore, anomalies in the Laplacian spectrum analogous to 
those of the HC network are observed for HMN matrices; in particu- 
lar, the lower edge of the HMN Laplacian spectrum has been recently 



shown to exhibit a continuous exponential Lifshitz tailp(/l) ~ e 

for with a ~ l'\ Taking the continuum limit of Eq.(4), we 

find p{t)~— / d?.p{?.)e^^'''-', which can be evaluated with the 

saddle-point method (see Methods), leading to 

p(f)~e-^, (5) 

i.e. anomalous stretched-exponential asymptotic behavior, in excel- 
lent agreement with computational results (see Fig. 4d). Therefore, 
hierarchical modular networks constitute a parsimonious and 
adequate model for reproducing all the complex synchronization 
phenomenology of the HC. 

A crucial role in the emergence of such behavior is played by 
disorder. One would be tempted to believe that all networks char- 
acterized by a finite spectral dimension could potentially give rise to 
this phenomenology. This is obviously not the case for a regular 
lattice, where the spectral gap is always well defined. A fractal lattice 
or an ordered tree, on the other hand, could exhibit a hierarchy of 
discrete low eigenvalues, whose multiplicities reflect system symmet- 
ries. The introduction of disorder, as in HMNs, is then necessary in 
order to transform such hierarchy of discrete levels into a continuous 
Lifshitz tail, leading eventually to the behavior predicted by Eq. (5). 

Discussion 

Simple models of synchronization dynamics exhibit an unexpectedly 
rich phenomenology when operating on top of empirical human 
brain networks. This complexity includes oscillatory behavior of 
the order parameter suggesting the existence of relatively isolated 
structural communities or moduli, that -as a matter of fact- can 
be identified by using standard community detection algorithms. 
Even more remarkably, oscillations in the level of internal coherence 



il l I T i Til I 1-1 i-l- I n I i-ffa 
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Figure 4 | Synchronization in hierarchical modular networks. Top panel: sketch of the HMN model. At hierarchical level 1, 2' basal fully connected 
blocks of size Mare linked pairwise into super-blocks by establishing a fixed number a: of random unweighted links between the elements of each (a = 2 in 
the Fig.). Newly formed blocks are then linked iteratively with the same a: up to level s, until the network becomes connected, (a), (b), (c) as in Fig. 3, but 
for a HMN with JV = 512, s = 5, and a = 4. Hierarchical levels are i = 1 ^ 5 in black, blue, green, magenta and red respectively (not all shown in a) for 
clarity), (d) Time relaxation of activity p for homogeneous characteristic frequencies to = 0, for logarithmically equally spaced values of k. Averages over 
Iff" realizations of HMNs with N = 4096 and s = 11. Inset: as in the main plot (d), but representing as a function of t"^ and confirming the predicted 
stretched exponential behavior, (e) Inverse tail-eigenvalues (as in Fig. 3) for a HMN as in e). 
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are also present within these moduh, suggesting the existence of a 
whole hierarchy of nested levels of organization, as also found in the 
recent literature relying on a variety of approaches' Aimed at 
unveiling this complex behavior we have introduced a family of 
hierarchical modular networks and studied them in order to assess 
what structural properties are required in order to reproduce the 
complex synchronization patterns observed in brain networks. 

In the absence of frequency dispersion, perfect coherence is 
achieved in synthetic hierarchical networks by following a bottom- 
up ordering dynamics in which progressively larger communities - 
with inherently different timescales- become coherent (see**^). 
However, this hierarchically nested synchronization process is con- 
strained and altered by structural bottlenecks -as carefully described 
here for the simpler two-block toy model- at all hierarchical levels. 
This structural complexity brings about anomalously- slow dynamics 
at very large timescales. Observe that the HC, in spite of being a 
coarse-grained mapping of a brain network, already shows strong 
signals of this ideal hierarchical architecture as reflected in its anom- 
alously slow synchronization dynamics as well as in the presence of 
non-degenerate eigenvalues in the lower edge of its Laplacian spec- 
trum, acting as a fingerprint of structural heterogeneity and com- 
plexity. We stress that such a complex phenomenology would be 
impossible to obtain in networks with stronger connectivity patterns 
(e.g. with the small world property) such as scale free-networks or 
high-degree random graphs. Even the generic presence of simple 
communities may not be sufficient to grant the emergence of frus- 
tration: the uniqueness of the human connectome, and of hierarch- 
ical modular networks in general, resides in the strong separation 
into distinct levels, which the synchronization dynamics is able to 
resolve only at well- separated values of the coupling k. 

On the other hand, in the presence of intrinsic frequency hetero- 
geneity, the described slow ordering process is further frustrated. 
Actually, for small values of the coupling constant k the system 
remains trapped into metastable and chimera-like states with traits 
of local coherence at different hierarchical levels. In this case, inter- 
moduli frequency barriers need to be overcome before weakly con- 
nected moduli achieve mutual coherence. This is clearly exemplified 
by the separation between distinct peaks in the chimera index x*" in 
Figs. 3-4, each one signaling the onset of an independent synchron- 
ization process at a given level (see Methods). The result is a complex 
synchronization landscape, which is especially rich and diverse in the 
intermediate regime put forward here. 

Including other realistic ingredients such as explicit phase frustra- 
tion** or time delays" *^ to our simplistic approach should only add 
complexity to the structural frustration effect reported here. It is also 
expected that more refined models -including neuro-realistic ingre- 
dients leading to collective oscillations- would generate similar 
results, but this remains to be explored in future works. 

Addition of noise to the Kuramoto dynamics would allow the 
system to escape from metastable states. Stochasticity can overcome 
the "potential barriers" between mutually incoherent moduli as well 
as re-introduce de-synchronization effects. These combined effects 
can make the system able to explore the nested hierarchy of 
attractors, allowing one to shed some light into the complex syn- 
chronization patterns in real brain networks. Actually, spontaneous 
dynamical fluctuations have been measured in the resting state of 
human brains^"; these are correlated across diverse segregated mod- 
uli and characterized by very slow fluctuations, of typical frequency 
<0.1 Hz, in close agreement with those found here (Fig. 3). 
Accordingly, it has been suggested that the brain is routinely explor- 
ing different states or attractors^' and that -in order to enhance 
spontaneous switching between attractors- brain networks should 
operate close to a critical point, allowing for large intrinsic fluctua- 
tions which on their turn entail attractor "surfing" and give access to 
highly varied functional configurations^'"""^ and, in particular, to 
maximal variability of phase synchrony"". 



The existence of multiple attractors and noise-induced surfing is 
largely facilitated in the broad intermediate regime first elucidated 
here, implying that a precise fine tuning to a critical point might not 
be required to guarantee functional avantages usually associated with 
criticality^^'^'''": the role usually played by a critical point is assumed 
by a broad intermediate region in hierarchically architectured com- 
plex systems'l Finally, let us remark that our results might also be of 
relevance for other hierarchically organized systems such as gene 
regulatory networks^*" for which coherent activations play a pivotal 
role™. 

Methods 

Numerical simulation of the Kuramoto model. The Kuramoto model is simulated 
by numerically integrating Eq. (1). Computations are carried out using both a 4th 
order Runge-Kutta method of fixed step size h — 10"' and an 8th order Dormand- 
Prince method with adaptive step size. Both methods lead to compatible results 
within precision limits. The robustness of the observation of a novel intermediate 
phase -between incoherent and coherent ones- is assessed by choosing different 
functional forms for the frequency distribution ^(w) (Lorentzian, Gaussian, uniform), 
and by implementing variations of Eq. 1 in which the matrix W is weight- 
normalized"'" for simulations of the HC and degree-normalized for simulations in 
HMNs. No qualitative change in the phenomenology is observed. 

Chimera index/*'' and hierarchical synchronization. In the main text, the chimera 
index y^^ is introduced as a measure of partial synchronization at the community level 
/. At any hierarchical level /, a hierarchical network can be divided into a set of 
communities. Following^^ is defmed as follows: (i) in the steady (oscillatory) state, 
and for each time t, local order parameters r|^' (t) for each community i are calculated 
and their variance across communities o'^h\(f) is stored; (ii) the chimera index is 
computed as the time average — ^cr^hi(f)^ ■ Having /''^ > 0 at a given hierarchical 

level / implies that local order is only partial as r-'' fluctuates, giving rise to a chimera- 
like state. On the other hand, — 0 means that each local order parameter at that 
level is rf ~ 1 , and local order has been attained. Figs. 3b-c and 4b-c show that at each 
/ (each color) a peak in the corresponding x"* marks the onset of the local 
synchronization processes: as soon as the peak vanishes upon increasing k, local order 
at that level is attained. The sequence of separated peaks in x'" for increasing values of / 
is the direct evidence of a hierarchical synchronization process. 

Lifshitz tail and stretched-exponential asymptotic behavior. In sparse HMNs, the 
lower end of the Laplacian spectrum is characterized by an exponential tail in the 
density of states p{l) ^ e" '^^ , known as Lifshitz taU'^. In graphs, Lifshitz tails signal 
the existence of non-trivial heterogeneous localized states governing the asymptotic 
synchronization dynamics at very large times t.ln the main text we have shown that in 
the absence of frequency heterogeneity, the f ^ cc behavior of the activity is given by 

p{t)^— / dAp{X}e^^^'-' . This expression can be evaluated by applying the saddle 

point method, yielding 

p(t)-yexp[-(l + fl)a-it;(2*:t)^]. (6) 

Substituting a ~ 1, as empirically found in HMNs^\ leads to Eq. (5), whose square 
root behavior is confirmed by simulations in Fig. 4d. 
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